function RF_filter_effect(wcdiff)

% RF_filter_effect shows the effects of the RF filter on the anti-aliasing filtering system
%
% RF_filter_effect(WC) plots the group delay and step response for the Analog Devices
% D68L8D filter block and the same block combined with the RF differential input lowpass
% filter. WC specifies the differential bandwidth of the RF filter in rad/s.

% Transfer function of the RF filter (differential response)
num_RF = wcdiff;
den_RF = [1 wcdiff];
sys_RF = tf(num_RF, den_RF);

% Transfer function of the predominant filter
[num_D68, den_D68] = D68L8D_tf(2*pi*250);
sys_D68 = tf(num_D68, den_D68);

% Transfer function of the cascaded systems
sys = series(sys_RF, sys_D68);

f = linspace(0, 250, 250);
w = 2*pi*f;
dw = w(2) - w(1);

[mag, phase] = bode(sys, w);
mag = squeeze(mag);
phase = squeeze(phase);
phase = myunwrap(phase);
delay = -diff(phase)/dw*pi/180;

[mag_D68, phase_D68] = bode(sys_D68, w);
mag_D68   = squeeze(mag_D68);
phase_D68 = squeeze(phase_D68);
phase_D68 = myunwrap(phase_D68);
delay_D68 = -diff(phase_D68)/dw*pi/180;

% Plot the group delays with an offset so that low-frequency ends overlap.
% Use the group delay of the combined systems as a reference since this corresponds
% to the delay one will see in practice.

RFtext = ['RF input-filter BW = ' num2str(wcdiff/2/pi) ' Hz'];
figure(1)
clf
f = f(1:end-1);
plot(f, delay, 'b', f, delay_D68 - delay_D68(1)+delay(1), 'r')
legend('RF cascaded', 'D68L8D module', 'Location', 'SouthWest')
title('Normalized Delay. RF cascaded response used as reference.')
text(5, 3.412e-3, RFtext)
ylabel('Time(s)')
xlabel('Freq (Hz)')

figure(2)
clf
step(sys, 'b', sys_D68, 'r', 10e-3)
legend('RF cascaded', 'D68L8D module', 'Location', 'SouthEast')
text(0.006, 0.1, RFtext)

% Find error at ~100 Hz
%[dummy, i] = min(abs(f - 100));
%[ f(i) phase(i)-phase_D68(i)]


function y = myunwrap(p_)
%  Unwraps degrees phases p_ by repeatedly looking for positive-going
%  discontinuities in p_ and subtracting 180-deg until the phase is monotonic.

p = p_;
y = p_;
loop = 1;
while loop
  % See whether there are any discontinuities in the phase
  dp = diff(p);
  i = find(dp > 0, 1);
  if ~isempty(i)
    y(i+1:end) = p(i+1:end) - 180;
    p = y;
  else
    loop = 0;
  end
end